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Dynamic Blind Signal Separation 

This invention relates to dynamic blind signal separation; that is to say blind signal 
separation which can be used inter alia in a non-stationary environment, and to a 
method, an apparatus and a computer program for implementing it. 

In recent years, blind signal separation (BSS) has emerged as an important field in signal 
processing. It is now successfully exploited in applications such as biomedicine, financial 
analysis, speech recognition and telecommunications. 

The objective of BSS is to recover or separate a number of statistically independent 
signals given only observations of mixtures of those sources detected using sensors, 
together with noise and unwanted interference. The term 'blind' is used to indicate that 
no prior information (e.g. waveform, frequency, bandwidth, location) is available 
concerning the individual sources or the manner in which they have been combined. 

■ 

Although the term 'blind* is often associated with these techniques there are several 
assumptions that need to be satisfied such that signal separation can be achieved. The 
mixing conditions, i.e. the combination process of the signals at the sensors, must be 
linear and time-invariant. The term time-invariant means that the spatial statistics 
(combinations) of the signals received at the sensor array do not change with time. For 
this invention, another assumption for the combination process is that the signals arrive 
at the same time at thQ sensors. This is referred to as the instantaneous mixing case. 

Many BSS methods have been developed for the problem when the combination 
(mixing) process is linear, instantaneous and time-invariant. These include: 

- The Blind Signal Separation (BLISS) algorithm, J.G. McWhirter, I.J. Clarke and G. 

Spence, "Multi-linear algebra for independent component analysis", SPIBs 44th 

Annual Meeting, The International Symposium on Optical Science, Engineering and 

Instrumentation, Denver USA, 18-23rd July, 1999. 
.... 

The Joint Approxirnate Diagonalisation of Eigenmatrices (JADE) algorithm with a 
description found in J.F. Cardoso and A. Souloumiac, 'Blind beamforming for non- 
Gaussian signals', lEE proc-F, Vol. 140, no. 6, pp. 362-370, Dec 1993. 



« 

"An Improved Cumulant Based Method for Independent Component Analysis", T. 
BlasGhke and L. Wiskott. Proc. Int. Conf. on Artificial Neural Networks, ICANN'02. 
2002. This reference relates to use of third order statistics and is an alternative to 
BLISS. It is a block based technique like BLISS and JADE. 

5 Published International Application No. WO 03/028550 A2 discloses a specim 
application of BLISS to monitoring fetal ECQ signals using abdominal eleptrodes applied 
externally to a mother. Electrode signals are filtered and subjected to a BSS 

r 

-technique based on Independent Component Analysis (ICA), I.J.Clarke, "Direct 
Exploitation of non-Gaussianity as a Discriminanr, EUSIPCO '98, Rhodes, 
10 Greece, 8-11 Septennber, 1998. 

When non-stationary conditions apply, the accuracy with which these existing BSS 
methods can detect and separate signals is restricted. This is because only a limited 
amount of data can be used reliably in the separation process, and statistically based 
estimates of the parameters may not be accurate. Additionally, processing of large data 
15 sets (e.g. over hours) is not possible. These limitations have led to the development o^. 
BSS methods for the non-stationary problem. Such methods include: 

ft 

- Particle filters for independent component analysis (ICA), R. M. Everson and S. J. 
Roberts, "'Particle filters for non-stationary ICA", Advances in Independent 
Components Analysis, Eds. M Gifolami, Sprinter, pp. 23-41 . 2000. 

20 - The EASI (Equivariant Adaptive Source Separation) algorithm. J. F. Cardoso and B. 
Laheld, "Equivariant adaptive source separation", IEEE Trans, on Signal Processing., 
vol, 44. no 12, pp. 3017-3030, Dec. 1996. 

- A Natural Gradient Algorithm (NGA) with Kalman filter approach, M. G. Jafari, H. W. 
Seah and J. A. Chambers, "A Combined Kalman Filter and Natural Gradient 

25 Algorithm Approach for Blind Separation of Binary Distributed Sources in Time- 
Varying Channels", IEEE conference on Acoustics, Speech and Signal Processing, 
Vol. 5, pp. 2769-2772, May 2001 . 

* 

- The Blind demixing algorithm, Z. Markowitz and H. Szu, "Blind Demixing Real-Time 
Algorithm of Piecewise Time Series". Internal Joint Conference on Neural Networks 

30 IJCNN'99, Vol. 2., pp. 1033-1037, 1999. 
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Particle filtering is a relatively new area of signal processing and it provides a robust 
* method for tracking signals. However, its computational cost renders its use impractical 
for real time BSS. 

* 

In contrast to this, the EASI and NGA (with a Kalman filter) algorithms are much more 
|||| 5 computationally efficient. These techniques- rely on stochastic and natural gradient 
approaches respectively. They tend to have slow convergence and in a non-stationary 
environment they are only suitable for tracking a slowly changing system. 

Alternatively, a more robust technique uses the 'block' algorithms (such as BLISS and 
JADE described earlier) but in a moving window approach. The term window is used for 
• 10 a section or 'block' of data. It is assumed that the data window advances in time to cover 
all recorded data as a series of sections and the signals are separated at each window 
location. These block algorithms process blocks or windows of data in isolation from one 
another. They are used in blind demlxing, and they assume that: 

- each window of data contains enough information about the signal mixing process 
15 and reliable estimates of the signals can be obtained; 

V 

- stationary mixing conditions apply locally, i.e. in each and neighbouring windows the 
mixing parameters will have a slow and smooth change. 

A major problem with using a moving window approach is called 'signal swapping'. This 
arises from unmixed signals being in different orders in successive windows without the 

20 orders being discernible: i.e. a first window might be processed to separate six signals 
numbered 1 to 6 in the order 123456, and processing of the immediately following 
window might yield these signals in a different unknown order such as 413562. The 
resulting observed signals therefore contain parts assembled perhaps randomly from 
different source signals. It may not be convenient, practical or possible to reorder the 

25 signals correctly to allow continuous monitoring of signals. 

The blind demixing algorithm overcomes this problem by assuming prior knowledge of 
signal statistics to reorder signals. However in many applications signal statistics- change 
from window to window and they may also be unknown. 

Moreover, as in particle filtering, the computational cost of a moving window technique 
30 can render its use impractical for real time BSS processing. 
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It is an object of the invention to provide an alternative method of signal separation. 

The present invention provides a method for dynamic blind signal separation including 
processing signals associated with data windows characterised in that the method 
includes: 

5 a) for pairs of successive data windows, each pair having a leading window anHf 
following window, using results obtained in connection with each leading vyindow to 
initialise data associated with the respective following window, and 
b) processing the initialised data to achieve signal separation. 

The invention provides the advantage that it is normally capable of processing data 
10 faster than the prior art of block algorithms, because initialisation avoids the need to 
carry out a full analysis of each data window (other than the first such) and it results in 
• computation being reduced. This is important in applications such as fetal heartbeat 
monitoring in particular, where processing in as near to real time as possible is desirable 
to enable clinical intervention in the case of an emergency. 



15 In a preferred embodiment, data orthogonality and data independence are bQ||| 
initialised, and they are initialised separately in respective stages. Data independent 
may be initialised by processing initialised data using independent component analysis 
(JCA) to apply small updates to the initialised data in a statistics procedure higher than 
second order to produce signal separation. 

* 

' 20 Data orthogonality may be initialised and the initialised data processed to update 
orthogonality using small updates to produce decorrelated data in a second order 
statistics procedure. The decorrelated data may be produced by a technique referred to 
as Jacobi and involving diagonalisation of a symmetric matrix by determining and 
applying rotations iteratively until off-diagonal elements of the matrix become 

25 substantially equal to zero. The method may include a second stage of initialisation using 
results obtained for each leading window to initialise independence of decorrelated data 
associated with the respective following window, this second stage using ICA to apply 
small rotation updates to initialised data in a higher than second order statistics 
procedure to produce signal independence and separation. The higher than second 

30 order statistics procedure may be at least one of a third order and a fourth order statistics 
procedure. 
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As will be described later in more detail, initialisation of both data orthogonality and data 
independence followed by small updates has the advantage of solving the problem of 
signal swapping. It also allows continuous processing of the signals. 

In a preferred embodiment, the invention is implemented in an acquisition phase in which 
5 signals are separated and desired signals are identified among the separated signals, 
and in a subsequent phase in which only desired signals are processed to separation. 
This has the advantage of greatly reducing the processing time required in situations 
where there are many signals to separate but only some (perhaps relatively few) of those 
signals are of interest, because apart from the acquisition phase unwanted signals are 
10 not processed significantly.- 

.In an alternative aspect, the present invention provides computer apparatus for dynamic 
blind signal separation programmed to process signals associated with data windows 
characterised in that the computer apparatus is also programmed to: 

a) process, pairs of successive data windows, each pair having a leading window and a 

15 following window, 

b) use results obtained in connection with each leading window to initialise data 
associated with the respective following window, and 

c) process the initialised data to achieve signal separation. 

Computer software for dynamic blind signal separation including processing signals 
20 associated with data windows characterised in that the software has instructions for 

4 

controlling computer apparatus to: • 

a) process pairs of successive data windows, each pair having a leading window and a 
following window, 

b) use results obtained in connection \yith each leading window to initialise data 
25 associated with the respective following window, and 

c) process the initialised data to achieve signal separation. 

The computer apparatus and computer software aspects of the invention may have 
equivalents of the optional and preferred features mentioned in connection with the 
method aspect. 

30 In order that the invention might be more fully understood, ah embodiment thereof will 
now be described, with reference to the accompanying diagrams, in which: 

Figure 1 is a schematic block diagram illustrating the method of the invention; 



Figure 2 shows signal scatter plots indicating stages in signal separation; 

Figure 3 is a block diagram illustrating the method of the invention; 

Figure 4 illustrates signal swapping in fetal EGG analysis; 

Figure 5 • is a vector diagram illustrating an update process employed in the inventio^^ 

5 Figure 6 is a scatter plot of two nearly independent signals; 

Figure 7 illustrates signal separation with the invention using maternal and fetal EGG 

data; 

4 

Figure 8 is a block diagram of the method of the invention adapted to target specific 

signals; 

10 Figure 9 shows EGG data prior, to signal separation and obtained with external 

- abdominal electrodes for a pregnancy involving twins; 

Figure 10 provides results obtained using the invention to separate signals from the 

Figure g data; 

« 

Figure 1 1 provides results obtained using the method of the invention adapted to target 
15 twin fetal EGG signals in the Figure 9 data; 

Figure 1 2 shows singleton fetal EGG data prior to signal separation; 

Figure 1 3 provides results obtained using the method of the invention adapted to target 

EGG signals in the Figure 12 data; 

Figure 14 shows pseudo three-dimensional plots of cross-correlated signals before and 
20 after correction for signal swapping; 

Figure 1 5 schematically shows electrode apparatus suitable for implementing fetal EGG 

monitoring; 

Figure 16 shows heartbeat/breathing data prior to signal separation and obtained, with 

sensors under a mattress on which a patient lay; and 

25 Figure 17 provides results obtained using the invention to separate signals from the 

Rgure 1 6 data. 



Referring to Figure 1. a dynamic blind signal separation (BSS) system 10 of tlie 
invention is arranged to monitor signals from two sources 11 and. 12. As illustrated, the 
system 10 has three sensors 14A. 14B and 14C (collectively 14) which develop electrical 
signals in response to reception of output signals 15 from both sources 11 and 12, 
although in actual applications of the invention more sensors may be used. One 
important application of the invention is monitoring fetal electrocardiogram (fECG) non- 
invasively: this uses sensors In the form of external electrodes applied to the abdomen of 
a pregnant mother. It is described in International Patent Application No. WO 03/028550 . 
(hereinafter "W550"), in which up to twelve signal electrodes are described defining 
respective signal channels in addition to ground and reference electrodes. It will be 
used as the basis for the present embodiment, which is concerned with the separation of 
a single fetal EGG signal for a singleton pregnancy, or multiple fetal signals for twins, 
triplets etc. 

The sensors 14 supply input signals to a data acquisition stage 16 of the system 10 in 
which the input signals are sampled in sets at a rate of approximately 500 sets of 
samples/second. The samples are then digitised and recorded. Each set of samples 
consists of a respective sample from each sensor 14, and the samples in each set are 
recorded substantially simultaneously. After recordal, these signals pass to a signal 
separation stage 17 where they are processed to yield separated or unmixed signals 
corresponding to source signals 1 5. The separated signals pass to an analysis stage 1 8, 
where they are analysed to distinguish fECG signal, maternal EGG signal and noise. 

The staicture of signals recorded In fEGG monitoring are complex, as the required fECG 
signal is subject to interference signals such as the much stronger maternal ECG signal 
and muscle noise. Typical prior art BSS methods are restricted to a relatively short 
section or time window of data over which prior art BSS assumptions hold good. There is 
however a well known and long felt want for signals to be separated and monitored over 
relatively longer periods of time, for which comparable prior art methods are not 
adequate. 

Data from the acquisition stage 16 is processed as a series of wiiat are referred to as 
data "windows" consisting of digital signals from each sensor 14, each signal consisting 
of 2,000-3,000 signal samples. Each window may partially overlap its preceding window, 
i.e. sets of samples may be common to both as described later. It is assumed that the 
data window advances through a succession of window locations in a bloci< of data, 
signals are separated for each window location and eventually the entire data block is 
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processed. A first data window is processed in the signal separation stage 17 by a block 
. mode method such as BLISS or JADE. Before describing this in more detail, a signal 
model for each window will be established and BSS used therewith. 

In fetal ECG analysis, the number of signal samples per set is normally taken to be eqy||^ 
5 to the number of signal electrodes, i.e. ignoring ground or reference electrodes. For BS^ 
an appropriate data model for data window number k having m digital signals from 
respective signal electrodes, each signal with n samples, is given by: 

X,=A,S, (1) 

where: is an (m x n) data matrix of window number k, and has rows which are 
10 respective signals (expressed as a set of m data samples); Sj, is an (m x n) data matrix 
of similar kind to Xj, but representing unmixed signals, i.e. signals which have been 
separated and corresponding to independent source signals 15; and A,, is an (m x m) 
mixing matrix which represents a transformation by which the source signals 15 are 
combined to form mixed signals received by electrodes and corresponding to X^,. 
15 mixing model in each window expressed by Equation (1) is assumed to be linear, 
instantaneous and time-invariant. Data in each window is assumed to have a zero mean 
(or indeed the data is zero-meaned). as this is convenient and does not affect the 
required result. 

The signal separation problem to be solved in BSS is to recover or unmix signals 
•20 represented by the data matrix S„: this carl be done by determining a matrix the 
inverse of Aj, , and premultiplying both sides of Equation (1) by it, i.e.: 

Ak'Xj,=A^^AkSk=Sk (2) 

Equation (2) for the kth data window can be written for the immediately following or 
(k+l)th window as: 



(3) 



To determine the inverse mixing matrix, a standard BSS algorithm can use two stages 



* • 
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The first stage, which only uses second order statistics, is to decorrelate and normalise 
the signals. The second stage, which normally relies on fourth order statistics, makes 
the signals (extracted from the first stage) independent. 

In order for these two processes to be carried out efficiency and effectively in each 
5 window, the tracking method has several important features. The first two include 
initialisation and update processes. 

The term initialisation means that information about the mixing matrix for the immediately 
• preceding data window is used to initialise signals associated with the current window. 
This initialisation has the effect of rendering initialised signals nearly separated, so that 
10 full recomputation is not required. To determine the mixing matrix of the current window, 
initialised signals will only need to be updated by a small amount. Moreover.- as will be 
described later in more detail, initialisation in combination with small updates can prevent 
the problem of signal swapping. 

It has been considered that the signals in the current window could be initialised using 

« 

15 A"' , such that a first approximation S'^^^ to S^^^ is provided by multiplying X^^^ by 



Aj^^ instead of A^^^ , i.e.: 



(4) 



The estimated S'^^, provided by Equation (4) will be closer to the desired value of S^^^ 
than X . However, in the two stage BSS approach (mentioned above) estimated 

20 Signals are normalised to have the same power. This has a consequence that in a first 
stage (decorrelation) signals cannot be updated. This will be shown later. Even if the 
normalisation process were to .be 'undone', updating the decorrelation process would 
then undo the effect of the previous stage. A feature of the example of the invention to be 
described is that it avoids this problem by updating the corresponding second and higher 

25 order components of the demixing matrix separately. 

In the present example, an initialisation process is implemented In two stages as.will be 
described in more detail later: in the first of these stages, the orthogonality of the signals 
is initialised and then the signals* orthogonality is updated using a decorrelation method. 
This makes the signals orthogonal. Then, in the second stage, the independence of the 
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signals is initialised and the signals' independence is updated using a higher order 
Statistics stage. This maizes the signals independent. 

In the first of the stages referred to above, one option Is to carry out a decorrelation and 
nornnaiising method on the (m x n) matrix X^, . This can be referred to as a singL^^ 

5 value decomposition (SVD) when the singular values are arranged in decreasing 
magnitude order (M. Moonen and B. De Moor, "SVD and Signal Processing, 111 
Algorithms, Architectures and Applications", ASIN: 0444821074, Elsevier 1995). 
However, ordering of signals in this manner gives difficulty, because signal powers can 
change from window to window resulting in signal reordering. Hence signal swapping 

10 can occur, i.e. the process generates apparent signals which are actually segments of 
different signals joined together. The terminology of decorrelation and normalisation is 
preferred here to SVD, as signals are not ordered according to their power. 

A second order decomposition of the matrix Xj^ can be expressed as: 

X, =u,2:,v; (5) n 

15 where XJ^is an (m x m) unitary matrix containing eigenvectors of X^X^; Vj, is an 
(n X n) unitary matrix containing eigenvectors of X^Xj^ ; superscript index T in IL^ or 

Vj" denotes a transpose of matrix Xj^ or ; and denotes an (m x n) diagonal 

matrix having singular values (square roots of eigenvalues) on its diagonal and zero 
values in all off-diagonal positions. These singular values are not necessarily arranged in 
20 order of decreasing magnitude and they indicate the relative contributions of associated 
eigenvectors to data to which both correspond. 

The columns of Uj, are respective orthonormal singular vectors: they are spatial vectors 

• indicating mixing of the source signals 1 5. The rows of Vj* are orthonormal singular 

vectors which are estimates of the source signals 15, but in general these estimates are 
25 not fully separated versions of source signals because they require a further relative 
rotation to become so. Decorrelation fails to separate the source signals fully because it 
is a second order decorrelation by which decomposed vectors are made mutually 
orthogonal, i.e. uncorrelated with one another. In many applications of the invention, the , 
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spatial characteristics (locations) of the signals will be similar, and so a solution that 
makes them dissimilar will not succeed in separating them. 

This failure to separate signals is illustrated in Figure 2, which is of the kind referred to as 
a "scatter plot". In this drawing a first signal is plotted on a vertical axis against a second 
5 signal on a horizontal axis in four plots (a) to (d). Points on the plots are indicated by 
squares and rectangles and represent values of the two signals at respective time 
instants (sample times). The plots (a) to (d) illustrate signals in various stages of mixing 
and unmixing. Plot (a) is a scatter plot of two source signals 15. which are assumed (for 
convenience of description) to be random uniformly distributed signals. 

10 Plot (b) is a scatter plot of two signals which are mixtures of the source signals, i.e. after 
mixing but before processing to unmix them. It shows that the mixing process 
corresponds to a complicated combination of stretching, shearing and rotation. Plot (c) is 
a scatter plot of two signals derived by decorrelation and nomialisatlon of the plot (b) 
signals. Such decorrelation and normalisation of signals is obtainable by application of a 

15 processing method which in statistical temns Is restricted to second order: this method 
provides estimated signals from which stretching and shearing have been removed. 
However, comparison of plots (a) and (c) shows that estimated signals in (c) are rotated 
relative to (a): i.e. the second order method has failed to remove a rotation. The 
estimated signals are therefore related to the source signals 15 by an unknown rotation 

20 matrix, and it is necessary to determine this matrix to rotate the estimated signals in plot 
(c) to obtain appropriately separated signals. Plot (d) shows two signals obtained by 
rotating those in plot (c). Processing which produces only a rotation does not produce 
shearing or stretching, and thus does not counteract results achieved by the second 
order method. The rotation is chosen to align edges of scatter plot (d) with co-ordinate 

25 axes (not shown). The correct rotation is determined as will now be described using 
higher order statistics, where "higher order" means higher than second order. 

In signal separation, detemrilnation of the unknown rotation matrix can satisfy a statistical 
independence condition. Mathematically, statistical independence of m signals Xi to 
X„ can be expressed as an ability to factorise the signals' joint probability function 
30 p(xi X ) into a product of the signals' joint probability functions p(xi) , i.e.: 
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m 



P(^i xj=YlM) (6) 



This process uses higher order statistics (HOS). The use of HOS to separate unknown 
and independent signals is often referred to as independent component analysis (IC^j^ 
and this is the second stage in the separation process in this example. 

5 From Equation (5), and only using data with non-zero singular values, Xi, can be 
expressed as Uj^ Sy, , where U,^ is an (m x m) matrix consisting of the first m 
columns of U^; is an (m x n)" matrix consisting of the first m rows of V,^ ; and Zy, is 
an (m x m) matrix containing the singular values of the kth window data matrix Xj, . From 
Figures 2(c) and 2(d), the signal vectors in are related to the source signals 15 by a 

10 rotation. Hence applying an'(m x m) unitary rotation matrix Rj, to will generate 
estimates of the unmixed or separated source signals 15: however, this will alter 
unacceptably the results obtained earlier by decorrelation and normalisation unless 

equal and opposite (m x m) counter-rotation matrix is also introduced at the same 
time such that R^Rj, =1, the identity matrix. Rl is applied to . This gives the 

15 following expression for X,^ 

X,=:U,Z^RXV,^ (7) 



In Equation (7), the product 2,, R^ is an estimate of the mixing matrix Aj^ and the 
product RfcV^ is a matrix having rows which are respective estimated signals. The 
estimated signals in R^V^ are orthonormal, i.e. orthogonal and normalised to have the 
20 same power, but have relative powers stored in the estimated mixing matrix . This 
normalisation aspect is inherent in the second order stage approach to the blind signal 
separation process. 



R^is calculated (as will be described later) by using an iterative sequence of pain^^ise 
rotations, each of which is designed to* maximise the statistical independence of a 
25 • respective pair of signals (taken from Vj^) as measured by a corresponding pairwise 
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cost function. In this maximisation process, the BLISS and JADE algorithms use fourth 
order statistics for the higher order statistics stage. In the present embodiment, the 
BLISS higher order statistics stage is used. 

For each window, the main objective is to calculate parameters contained in Equation 
5' (7). If each window were to be processed in isolation as in the prior art. then results for 
each window would be determined separately: the computational cost of using k 
windows would then be k times that for a single window. This is impractical if the number 
of windows is large, e.g. > 20 windows with 2,500 samples per signal per window. 
However, as previously described it has been found that it Is possible to reduce the 
10 processing required by using results obtained for each window In initialisation of the 
separation process for a subsequent window. This initialisation converts data into an 
intermediate form which is closer to the desired result and reduces processing. It results 
in signals which are closer to being separated and a full recomputation of signal statistics 
Is not required. The mixing matrix need not be stationary to permit this, but should 
15 be changing relatively slowly with time. 

Having discussed the BSS problem, it can now be shown that initialising signals with A,, 
is not feasible. For example, rewriting Equation (7) by substituting index (k+1) for k, 
multiplying both sides by (u^Sj^R^^j . the inverse of the estimated mixing matrix 
A for the kth window, and replacing inverses of matrices by transposes::- 

20 (tJ, K T = (R. U: = (Ru )K^ S... RLxR.« Vj« ) (8) 

In Equation (8), (u,« RLiR^^Vv*,) is unknown and has to be estimated from 
the data. U^and R^are unitary matrices, and therefore U;^ =tJl, XJlV^ =1 and 
Uj.tJk = I . where I denotes the identity matrix. If A,, is slowly changing, then 
UkP^^ - 1 . RkRLi - 1 ^<^« " ^ • '-'^'"9 ^^^^^ approximations:- 

25 (Uk£.R;[)"x,«»R,^,V,\. (9) 
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On the right of Approximation (9), the product I^+iX^^i ^ matrix having rows which are 
respective estimated signals: these signals are initialised and close to orthonormal, i.e. to 
orthogonality and equality of signal powers, and can be processed by moving window 
ICA using a decorrelation method. However, they were obtained using a .decorrelation 
5 and normalisation method in Equation (5), and this discriminates between unnormali^^ 
signals on the basis of signal power. Such discrimination makes Equation (5)*s 
decorrelation process much more difficult. A data window initialisation process is 
preferred (as will be described) which initialises the orthogonality of the incoming signals, 
i.e. those in the (k+1)th data window, but which avoids normalisation of these signals. 

io The present example is. implemented in two stages, each- involving a respective 
initialisation process: in the first of these stages, the orthogonality of the signals is 
initialised by applying the kth window matrix . to them, and then their orthogonality is 

updated to reflect data in the (k+1)th window using a decorrelation method. This makes, 
the signals orthogonal. Then, in the second stage, the independence' of the signals is 

15 initialised by applying the kth window rotation matrix Rj, to them, and then their 
independence is updated using higher order statistics. This makes the sign 

* 

independent. 

Figure 3 is a block diagram of signal separation in the present example. It summarises 
this example, and to relate it to theoretical analysis, reference will be made to equations 
20 to be provided later. Processing of the kth data window yields initialisation information in 
the form of matrices Uj, and Rj, . The matrix Uj, is input at 21 and used at 22 to 

initialise orthogonality of data in an immediately following (k+lth) window using Equation 
(16). Initialised data are decorrelated at 23 using Equation (18) with small update angles, 
i.e. "smair means more than one possibility is available and the smallest angle is 

25 chosen: step 23 uses a technique referred to as "Jacobi" and described in Adaptive Filter 
Theory (3rd Edition), S. Haykln. 1996. The expression "update" is employed because it 
reflects the fact that step 23 adjusts the results obtained at step 22 so that they are 
updated from partial reliance on the kth data window to become more fully reliant on the 
subsequent (k+1)th data window. Jacobi involves diagonalising a covariance matrix of 

30 the data matrix X^^^ for the (k+1)th window, the covariance matrix being a symmetric 

matrix. Steps 22 and 23 collectively form a second order stage of processing in statistical 
terms producing orthononmal (orthogonal and normalised) signals. 
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The matrix is input at 24 and applied to the orthononnal signals at to initialise them 

using Equation (26). The initialised orthonormal signals then undergo separation by 
rotation at 27 by ICA with small angle updates using higher order statistics and Equation 
(29). Here again, "small" means more than one possible updating angle is available and 
5 the smaller/smallest angle is chosen. Steps 25 and 27 are collectively a higher order 
stage of processing in statistical terms producing separated signals at an output 28. 

In many decorrelation routines, signals are made orthonormal by using an iterative 
sequence of painwise rotations. Each of these rotations is designed to maximise the 
orthogonality of a given pair of signals as measured by a corresponding painwise (second 
LO order) cost function. In the Jacobi method, this can be accomplished by reducing a 
symmetric matrix to diagonal form using Given's rotation as described in the Haykin 
reference mentioned above. For example, in a two signal case the objective of using the 
Jacobi method is to diagonallse a (2 x 2) covariance matrix X^^ of two sianals Xi and Xa 
expressed as (n x 1) vectors. This matrix is given by Haykin as: 

(10) 

This can be implemented by iteratively forming a (2 x 2) matrix B by. rotating the Equation 
(1 0) covariance matrix using a rotation matrix Q to force off-diagonal elements of B to 
zero, i.e.: 
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B = 



12 



t>21 ^22 



= Q 



T T 

Xj Xj Xj^ X2 

T T 

XjXj X2X2 



Q = 



COS Q - sin 0 
sin @ cos G 



T 

X^ Xj 



T 

XjXj 



T T 
JC2^l X2X2 



cos © sin 0 
- sin © cos © 



(11) 
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30 



At each iteration of Given's rotation, the rotation matrix Q is determined which minimises 
off-diagonal elements of the symmetric matrix • When the matrix becomes fully 
diagonal, i.e. when ail its off-diagonal matrix elements have become zero, the signals it 
represents have become mutually orthogonal. This can be achieved by simply setting 
b to zero in Equation (11), and determining a value of 0 that satisfies this. B is a 
symmetric matrix, so bij= bj^ and it is only to determine 0 for h^^= 0. By forcing h^^ 



t 

f 
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and b^i to be equal to 0, and Xa are nnade orthogonal and hence B becomes a 

diagonal matrix. Thus, the diagonalisation of B can be a necessary condition for making 
vectors orthogonal. 



Using Equation (7), a covariance matrix Xj^X^ is formed for Xj^ , i.e.: 




Because Rj^ is a unitary matrix, and because estimates of signals contained in the rows 
of Vk are orthogonal to one another: 



(13) 



Then X^X]^ = U, £ =Ukaktj;[ (14) 




10 where is an (m x m) diagonal matrix whose diagonal elements are eigenvalues of 



Xj^. 



Writing (k+1) for k in Equation (14), for the (k+1)th window the covariance matrix 



Xt.,Xt., is expressed as : 



15 In order to avoid computing decorrelation afresh, the covariance matrix ^m^m-^^h-^i 

initialised by premultiplying by and postmultiplying by , which were obtained for 
the immediately preceding or kth window, i.e.: 

r 

. Equation (16) represents tJj, initialising orthogonality for information or signals in the 
20 data matrix Xj,^! using results obtained for the kth data window. This initialisation is at 
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22 In Figure 3: it improves the ortliogonaiity of the vectors in the covariahce matrix. It 
also makes the matrix U^(X,^,X;[,i)U, on the left hand side of Equation (16) close to 
being diagonal, i.e. off-diagonal elements of this matrix become small. As previously 
mentioned, the covariance matrix being diagonal is a condition for making vectors in that 
5 matrix orthogonal. Initialisation at 22 is followed by a further improvement in orthogonality 
brought about using information from the (k+1)th data window as described below. 

Now in a slowly changing system, U][U^„ = I and U^«Uk = I. and so Equation (16) 
can be rewritten as: 

u^(u,cA.iUL)Ui: = a,« (17) 

10 For an (m x m) symmetric matrix, and without initialisation of data, a Jacobi 
diagonallsation method is iterative with of order m^ operations required for each sweep 
consisting of a single update of ail the relevant signals represented as rows of the matrix 
X X3,. Many sweeps may be required. However, the matrix (UkX,j+iXk+iUk), 
initialised at 22 in accordance with the invention and given by Equations (16) and (17), Is 

15 already close to representing a diagonal matrix: the Jacobi diagonalisation method 
therefore requires only a few iterations to converge. Thus, the computational cost Is 
reduced in accordance with the invention by using an transformation to initialise signal 
orthogonality in the (k+1)th window, i.e. the current window. 

The next step Is to complete diagonalisation for the (k+1 )th window using in this example 
20 the Jacobi method at 23. It employs matrices and which are eigenvectors 
estimated when applying Jacobi to the matrix (UkXk^.iXj+iUk) which has been , 
initialised as regards orthogonality. and tj^ are rotation matrices implementing 
small rotations that diagonalise the initialised matrix (UkXk+iXk+itJk). i e- they should 
make vectors in this matrix orthogonal by updating as follows, i.e. : 

25 tJl(u][x,«xLA)u, = a,,, (18) 

In Equation (18), U, is computed using Jacobi- diagonalisation of the orthogonality 
initialised covariance matrix in Equation (1 6): a unitary matrix is required to diagonalise 
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this matrix X^^,Xj[^j . As before, a transform in the form of a rotation Implemented by a 
matrix Q is sought such that (x,,^XJ,Jq is diagonal, where Q =U,U . This 
uses the fact that U,, initialises data and is based on small updates. They can be 
used in the form of their product as they are both are unitary and their product is t 
also unitary. Estimates of eigenvectors t^^, for the (k+1)th window are expressed by: 



and the estimated signals ( V.^^J are defined in Equation (24). This all occurs in stages 
22 and 23.' 

For the higher order stage 25/27. the signals obtained from Equation (24) are initialised 
in Equation (26) using . This occurs in stage 25. shown in Figure 3. This process 
initialises fourth order information in the signals. To make the signals more fully 
independent or separated, they are updated using Equation (29). This occurs in stage 



27. 



Writing k+1 for k in Equation (5): 



X =U y 

"^k+l k+1 -^k+l ^ k+1 



(20) 



Premultiplying Equation (20) by ■ 



k+l-^k+l ^k+l~-^k+l '^k+1 (21) 



Premultiplying Equation (21) by : 



k+i (22) 



^k+i ^k+iX^^j - Efcli V^^, = VJ 



where is a diagonal matrix having elements which are the square roots of the 
eigenvalues contained in the matrix a^^, , i.e. 



t 
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(23) 



Orthonormal signal vectors contained in the matrix V^j (for use in the higher order 

0 

Statistics stage 25/27 of the method 20) can be expressed by eliminating Sj^+i from 
Equation (22) using Equation (23), i.e.: 



^ k+1 "^k+l ^ k+1 k+1 



(24). 



The vectors or signals Vj^^j are decorrelated signals output from the second order stage 

22/23 in Figure 3. As described previously, they are related to unmixed signals 
(corresponding to independent source signals) by an unknown rotation matrix, as 
described with reference to Figures 2(c) and 2(d). The higher order statistics stage 25/27 
10 is used to estimate this rotation matrix. To reduce computation once more, information 
from a previous window is used to initialise the separation process for a current window. 

For the (k+1 )th window, m orthonormal signal vectors each with non-zero eigenvalues 
and zero a,ean have bean obtained as V^, from *e Jacobl decorrelafcn stage 23. 

They are now expressed as a product of an (m x m) unitary rotation matrix Rj^j and 

15 signal estimates Sj^^.^ which are unmixed versions of the source signals 15.. 



Vfc+i =Rk+iSk+r (25) 



The estimates Sj^+j are related to the source signals 15 by an arbitrary permutation to 
reflect possible signal swapping and a scale factor. It is necessary to estimate a rotation 
matrix R^+i which is the transpose and inverse of R^+j in order to multiply both sides of 

20 Equation (25) by it: this eliminates Rl^^ from the right hand side of Equation (25) and 

' consequently recovers the required estimates S^^^. If the mixing matrix is changing 
slowly, then by using information from the preceding or kth window: 

R^V^^+i =RkRk+i^k+, (26) 



. 20 

Equation (26) expresses decorrelated signals output from step 23 having their 

independence initialised at 25 by multiplication by Rj, at 25 obtained from the kth data 

Window. This initialises fourth order information in the decorrelated signals. To improve 
signal independence and separate the signals, they are updated in step 27 on the basis 
5 of the )th data window as follows: 

R^RLi ==1 for slow change, and so: 





R.V,^«=S,,, (27) 

> 

The higher order statistics stage 25/27 in a block based algorithm such as BLISS is 
implemented in accordance with the invention using the initialised matrix R^Vj^^.! of 
10 signals from Equation (27). This stage is iterative with approximately the following 
number of operations per sweep, a sweep being a single update of all signal pairs in the 

matrix R^Vj+i . 

No. of operations per sweep « 20m(m-1)n/2 (28) 

. where the term m(m-1)/2 denotes the number of pairwise operations and n denotes the 
15 number of signal samples in the (k+1)th window or elements in the associated data 

matrix Xj^^i . Since the signals defined in Equation (27) are close to being separated. 

only a few iterations of sweeps in the higher order statistics stage 25/27 will be required 

» 

for convergence of the signals to separation. Separation is judged to have occurred 
when when angles obtained from the BLISS estimation process are small. These angles 
20 are estimated from Equations (35) and (36) to be described later. Initialisation in 
accordance with the invention reduces the number of computation steps required to 
unmix signals. Signal independence is improved by updating in step 25, and unmixed or 

separated estimates S^+i of original signal sources are given by: 

S,,,=R,(rA"i) (29) 

25 where R^ is a unitary rotation matrix implementing a small rotation or update. R^ is 
computed using the initialised vectors (RkVj^+i) in the higher order stage 25/27. Similarly 
to R^ , it is calculated by using an iterative sequence of painA^ise rotations, each of which 
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is designed to maximise the statistical independence of a respective pair of signals 
(taken from V^+j , as measured by a corresponding- pairwise cost function. A unitary 

rotation matrix Rj^^^ for the (l<+1)th window is required to rotate signal vectors extracted 

* 

by decorrelation as described earlier. This uses the fact that Rj, initialises data and R^ 
5 is a small update to Rj, : they can be expressed as a product as they are both are 
unitary and their product is also unitary. The rotation matrix Rj^^^ is therefore expressed 
as a product of an initialisation matrix and an update matrix, i.e.: 

R.« = RzR. (30) 

The initialisation and update processes described with reference to Equations (19) to 
10 (30) can reduce the computation load associated with separating signals. Additionally, 
they can be used to prevent a known problem referred to as signal swapping as will be 
described later. 

To summarise, Equation (16) represents initialisation by Uj, of orthogonalisation of the 
information (vectors or signals) in the data matrix X^^^. This corresponds to step 22 of 
15 Figure 3. To make the signals more fully orthogonal, they are updated in step 23 using 
Equation (18). Updated eigenvectors (tJ^+i) are estimated using Equation (19) and 

'estimated signals {^^^^ ) are defined in Equation (24). For the higher order stage 25/27, 
at 25 the signals obtained using Equation (24) are initialised by Rj, using Equation (26). 
Step 25 initialises fourth order information in signals being processed using it. To make 
20 these signals more fully independent or separated, in step 27 they are updated using 
Equation (29) to more fully reflect the current (k+1)th window. 

Attention will now turn to the problem of signal swapping, which is illustrated in Figure 4 
for fetal EGG analysis. Data windows were used which were 4 seconds/2,000 samples 

« 

long and 8 samples wide: i.e. there were eight signals each of 2,000 samples per 

« 

25 window. A window overlap of 50% was used: i.e. the second half of each window was 
the first half of the immediately following window. The windows were processed in 
isolation and the sensors 14 consisted of eight signal electrodes together with ground 
and reference electrodes. This drawing shows eight signal traces 1 to 8 (marked, on left) 
of amplitude against time: they are each intended to represent an individual separated 

30 signal, but instead each contains contributions from multiple signals. For example, trace 
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3 contains a first region 30 with peak periodicity appropriate for fECG, a second region 
32 with peak periodicity appropriate for maternal ECG and third regions 34 of noise. 

The phenomena of signal swapping can potentially occur in the second order stage 
22/23 (Equations (1 9) to (22)) and/or the higher order stage 25/27. It is counteracted 
applying only small -^4 < 8 < +7i/4 rotation updates in the initialisation process previoi 
described. For the purposes of the discussion to follow, x and y are vectors and each is 

a respective row of Vj^ representing all samples of a signal from a respective sensor 
extending the full length of a data window. After x and y are initialised in the second 
order stage 22/23, they are close to being orthogonal to one another. This is shown in 
Figure 5, where initialised vectors x and y (solid arrows) are close to 90 degrees apart 

as indicated by orthogonal co-ordinate axes indicated by dotted lines 40. 



To make the vectors x and y more accurately orthogonal, a relative rotation angle 6 is 

introduced between them. This rotation is implemented by a decorrelation method such 
as Jacobi diagonalisation of a symmetric matrix as described earlier. For signals which 
15 are already close to orthogonal, two solutions are obtained by this method, i.e. 0 0 a 
0 « 90 degrees. To avoid signal swapping the smaller value of Q^O should be chosen. 
If the signals are indeed rotated relative to one another by the smaller value of 0, this 
corresponds to a Givens rotation matrix given by: 




' cos 0 


sinO 






0' 


-sin0 


COS0 




0 


1 



(31) 



20 When 0 = 0 then a pair of initialised (nearly orthogonal) vectors x and y can be 
updated using: 



r -Ti 
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.T 

.y 




0 1 


T 



(32) 



In this case, the updated signals are in the same order and signal swapping, has been 
avoided. This is in direct contrast to the case of 0==9O, for which a Givens rotation 
25 matrix for rotating one vector relative to another is: 
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- sin 6 cos 6 




-1 0_ 



(33) 



Updating a pair of initialised vectors x and y using Equation (33) gives: 
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(34) 



■ i.n this case, the updated signals have swapped as indicated by the exchange of 
5 positions of x"^ and -y"^ and there also occurs an inversion indicated by the negative 

sign of - y ^. 

In the Jacobi decorrelation method, the smaller of the two angles is therefore chosen to 
achieve convergence to a solution. Thus, signal swapping is prevented in accordance 
with the invention when using an initialisation process for the second order stage and 
10 . Jacobi with smaller angle updates. 

Signal swapping can also potentially occur in the higher order statistics stage, as again 
there will be more than one solution or rotation angle to choose from. After is applied 
to data using Equation (27). estimated signal vectors will be related to corresponding 
source" signals by a small rotation. In the BLISS or higher order statistics stage, since 

15 fourth order statistics are used., there are now four angles that can be chosen to achieve 
signal separation. This is illustrated in Figure 6. which shows a scatter plot 50 of two 
nearly separated signals rotated through a small angle relative to Cartesian co-ordinate 
axes 52. If these signals were properly separated, the lower and left hand edges 54 and 
56 of the scatter plot 50 would be aligned with respective axes 52. When signals are 

20 nearly separated in this way. the scatter plot 50 can be rotated by any one of four angles 
(^ = 0,90, 180 or 270 degrees) to align it with axes 52. Here again, in accordance with 
the invention, it has been found that choosing the smallest of the four angles (^ = 0 
degrees) avoids signal swapping. 

The estimation process in BLISS does not automatically choose the smaller, angle, yet 
25 signal swapping does not occur. BLISS attempts to estimate the painwise rotation angle 
between signal axes and axes of strongest two-fold or four-fold symmetry in a two- 
dimensional scatter plot. Four-fold symmetry is best suited to the case Where both 
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Signals have a distribution with l<urtosis (fourth order statistics) of the same sign. Two- 
fold symmetry is most appropriate when the signals have kurtosis of opposite sign. 

In the case of four-fold symmetry, for a pair of orthonormal (zero-mean) n-element 
vectors x and y representing signals in digital form, the pairwise rotation angle 9^ 
estimated by the formula: 



49 = arg ument 



v. 



(t) cos4(t)(t) + j^r'^ (t) sin 4<j)(t) 

J 



(35) 



where :x(t). y(t) = x. y respectively at time t, (j) = tan'' y(t)/x(t) , =x^(t) + y2(t). 
x(t) = r(t) cos (j)(t) , y(t) = r(t)sin ^(t) and = -1 . 

In the case of two-fold . symmetry, for a pair of orthonormal (zero-mean) n-element 
10 vectors x and y . the palnwise rotation angle is estimated by the formula: 



15 



20 



25 



49 = arg timeiit 



2,r (t)cos2<t)(t) + j2r''(t)sin2<|)(t) 
t=i t=i J 



(36) 




In general a combination of the two estimators may be used .to good effect. 

In the BLISS higher order stage the scatter plot is actually rotated by 46. For nearly 
independent signals this means that the scatter plot is simply rotated by full cycles (i.e. 
by (0.90.180.270) degrees multiplied by 4). It then follows that the estimated 40 solution 
will be close to zero. For example, the higher order stage in the BLISS algorithm uses a 
cost function that has sine terms in the numerator and cosine terms In the denominator. 
The sine of 6. when 9 is approximately given by 4x(0.90.1 80.270). Is close to zei-o and 
the estimated angle in BLISS will also be close to zero. Thus, the same principles that 
applied for the initialised signals in Equations (31) - (36). will also apply in this case. 

As has been said. BLISS is based on fourth order statistics. It is also possible to use a 
technique based on third order statistics instead of BLISS, as described in "An Improved 
Cumulant Based Method for Independent Component Analysis". T. Blaschke and L. 
WIskott. Proc. Int. Conf. on Artificial Neural Networks, ICANN'02, 2002. It is also 
possible to use a combination of third and fourth order statistics. 
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The invention as described above was used for fetal EGG analysis using the data shown 
analysed in Figure 4 but with signal swapping. Data windows were used which were 4 

, • 

seconds/2,000 samples long and 8 samples wide with a window overlap of 50%: i.e. the 
second half of each window was the first half of the immediately following window. The 

5 sensors 14 consisted of eight signal electrodes together with ground and reference 
electrodes to yield eight signals. The windows were processed as indicated in Equations 
(19)-(30). Results are shown in Figure 7. This drawing shows eight traces 61 to 68 
(marked on left) of signal amplitude against time: they can clearly be seen to represent 
well-separated signals because signal characteristics do not change along each trace. 

10 For example, traces 61 and 65 both show ~ 140 beats/minute appropriate for fECG; 
traces 62. 67 and 68 show ~ 87 beats/minute appropriate for maternal EGG. A noisy low 
frequency signal appears well separated on trace 64. and noise Is separated out at 63 
and 66. This demonstrates the ability of the invention to prevent signal swapping. The 
reason for there being e.g. more than one maternal EGG signal is that the heart is not a 

15 point source and there may be more than one EGG signal from the same extended 
source. 

The initialisation and update processes also allow desired signals to be targeted alone. 
• This technique avoids redundant computation incurred by tracking unwanted sig^nals. It 
exploits the fact that signals of interest (e.g. maternal EGG and fetal EGG) can be 
20 classified in an acquisition stage. For fetal EGG analysis, this classification process could 
either depend on the choice of the user (clinician or doctor) or on a pattern classification 
technique. Pattern recognition methods for EGG signals are described by J. Pan and W. 
J. Tompkins in "A Real-Time QRS Detection Algorithm". IEEE Transactions on 
Biomedical Engineering 32(3): 230-236, 1 985. 

25 Once the desired signals are identified, this information can be used to focus or target 
processing on them. After classifying signals bf interest In the acquisition stage 16. as 
compared to the invention without such targeting, the targeted technique only differs in 
restricting updating to desired signals in the signal separation stage 17. 

Figure 8 Is a block diagram of signal separation 80 with desired signal targeting. 
30 Initialisation Information consisting of the matrix Uj, from an immediately preceding (kth) 
data window Is input at 71 , and at 72 is applied to data in the (k+1)th data window 
to initialise its orthogonality. Desired signals are.then targeted at 73 as will be described 
in more detail later, and only these desired signals (not others) are decorrelated at 74 
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with small update angles using Jacobi. Steps 72 to 74 are collectively and in statistical 
terms a second order stage of processing producing ortHonormal signals. The kth 
window' rotation matrix is then input at 75 and applied to the orthonormal signals at 

76 to initialise their independence. Desired signals are then targeted once more at 77 
5 The desired signals undergo separation at 78 by ICA with small angle updates us 
higher order statistics. Steps 76 to 78 are collectively a higher order stage of processing 
in statistical terms producing separated signals at an output 79. 

Targeting of desired signals at 73 and 77 in Rgure 8 appears in both the second and 
higher order stages. This Is not essential, one or the other may be used instead of both. 
10 Most benefit will be gained by targeting desired signals at 77 in the higher order stage. 
This higher order targeting approach will again use a process of initialising second order 
information, updating with Jacobi using Equation (1 6) and then initialising higher order 
data using Equation (27). 

Before describing the targeting method in more detail, it is re-emphasised that a higher 



rotations. 



age 76 to 78 can be implemented using an iterative sequence ot pairwisa|^ 

I. Each of these rotations is designed to maximise the statistical independence' 

of two signals forming pair of signals, statistical independence being measured by an 

associated pairwise cost function. With m signals to \^ , .there are "^€2 or m(m-1)/2 
signal pairs. After initialisation of these signals to form orthonormal vectors ^\ to v'^ 
20 with zero mean at 76, the total number of vector pairs are expressed as: 

vX v;v; v;v; ••• vy^ 



. : (37) 



/ 

m 



In (37), it is assumed that the orthonormal vectors y\ to are initialised using higher 

order information obtained using an immediately preceding data window as described 
• previously. All of these pairs would have to be updated in an iterative sequence in the 
25 absence of targeting, which means that m(m-1)/2 pairs would need updating. In contrast 
to this, the targeting approach requires only to update pairs that include the desired 
vectors. This is possible as: 
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the initialisation and small update processes allow the vectors to maintain their order 
in the separation process (for each data window their locations in the triangle will not 
change); and 

■ by identifying the desired signals (maternal and fetal ECG in this example) in a signal 
acquisition stage, these signals can be denoted as vectors and in updating. 
Pairs involving only v[ and individually and collectively can then be targeted. 

If pairs Involving and Vj are targeted and no otliers. only the first two rows of (37) 
need updating. This involves only (m-1) + (m-2) signal pairs.' as opposed to m(m-1)/2 
pairs for a full update. In this targeted* process, only the maternal and fetal signals will 
be separated. Separation of other signals is not accomplished, but these are unwanted in 

any event. 

As described above, the targeting principle can also be applied to the second order stage 
72-74. After an acquisition stage and by initialising second order information as in 
Equation (11) to produce an initialised covarlance matrix, rows of this matrix that 
correspond to the desired signals will be updated, but not other rows. 

Referring now to Figure 9, there is shown ECG data obtained from a mother pregnant 
with twins. This data was obtained using 1 6 external abdominal signal electrodes with 
respective signal channels together with ground and reference electrodes as described 
in W550, except that signal separation was not used. Data was gathered over a 30 
second time interval. There are 1 6 signal traces 91 to 1 06 ( marked on left), in most of 
which (92, 93. 95-1 06) the strong maternal ECG signal can be seen as spikes such as 
107 spaced by about 0.8 seconds (-72 beats/second). Fetal ECG signals from the twins 
are not discernible over noise. 

Figure 1 0 shows results obtained using the technique of the invention to separate signals 
shown in Figure 9, but without selecting out desired signals by targeting as described 
above. To process these signals, data windows were used which were 6 seconds/3000 
samples in length and 1 6 samples in width because of the 1 6 signal channels. Each 
window overlapped the respective immediately preceding window by 60% (3.6 seconds). 
The drawing shows signal traces 11 1 to 126 ( marked on left), where traces 1 12 and 125 
are the maternal ECG. Traces 1 15 and 1 18 are the fECG signal for one twin, and traces 
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121 and 123 are the fECG signal for the other twin. The traces show that the signals are 
well separated and continuous, i.e. signal swapping has been prevented. 

Application of the invention to produce the results shown in Figure 10 involves the 
computational burden of tracking 16 signals. A signal detection procedure, such as 

5 use of singular values in the acquisition stage, could have been used. Singular val 
provide an indication of the power of the signals in the data. However, as a fetal EGG 
signal is often weak this detection method is not accurate. Nevertheless, on separating 
signals in an acquisition stage they can be Identified (using pattern matching) and the 
method of the invention can then be targeted to focus on desired signals. Pattern 

10 recognition techniques can be used to look for repeating shapes: e.g. a heart beat signal 
has a central section with a characteristic sihape referred to in cardiology as a QRS 
section. The QRS of the maternal EGG should be much stronger and repeat less often 
than the QRS section of the fetal signal) 

Figure 1 1 shows the results of applying the method of the invention in targeted form to 
15 fetal EGG signals from respective twins using the same d^ta as that used in relation to 
Figure 10. Two signals 140 and 142 are shown, each from a respective twin, and 
again continuous, i.e. signal swapping has not occurred. The targeted method used the 
same acquisition stage 16 as the untargeted method. Signals 115 and 121 in Figure 10 
were identified by inspection as fetal EGG signals of different twins and these were 
20 tracked using the method of the invention In targeted form to give the results shown in 
Figure 1 1 . In the targeted approach, a first block of data is used in an acquisition stage, 
and all available signals are processed to separation so that desired signals may be 
selected from them. Desired signals may be identified either by: 

■ the user by inspection, or 

25 ■ one of a number of pattern recognition techniques, which are widely available in 
scientific literature. 

After the acquisition stage, only desired signals are processed to full separation. Once 
the order of the separated signals has been determined in the acquisition stage, as 
• shown in the triangle of signals (37), the signals can be re-ordered such that the desired 
30 signals occupy the top rows of the triangle: i.e. signals 115 and 121 in Figure 10 would 
become the first and second signals in that example. The initialisation and update 
processes then maintain this order, i.e. the signals can be targeted because their 
positions in the update process are known. 
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In a second example, the targeted naethod of the invention was applied to a singleton 
fECG recording, i.e. an EGG of a single fetus. The recording employed twelve signal 
electrodes and associated ground and reference electrodes and lasted for three minutes. 
In Figure 12, a thirty second Interval of the recording is shown of twelve signals (a) to (I) 
associated with respective electrodes. To process this recording of a data set, data 
windows of length 6 seconds with zero window overlap (i.e. contiguous windows) were 
used. The results are shown in Figure 13 as three plots against time for the first (0-1), 
second (1-2) and third (2-3) minutes of the recording respectively, and like signals are 
like-referenced. The plots show signals 150 and 152, which correspond to the fetal and 
maternal signal respectively. These signals are continuous, I.e. they have not swapped. 

Another method for dealing with -signal swapping is to use overlapping windows and a 
cross-correlation of overlapping estimated signals which result. For example, let 
overlapping sections of estimated signals in windows k and k+1 be f samples long and 
be represented by (m x f) matrices and S^^^ respectively. The rows of these 
matrices contain the estimated signals and the overlapping sections will contain similar 
signals. The niatrix S^is related to S^+^ by: 



S^,,«DPS^ . (38) 



where D denotes an (m x m) matrix, whose diagonal elements are 1 or -1 (providing for 
signal inversions) and P denotes an (m x m) permutation matrix (providing for signal 
swapping). These conditions arise in adjacent windows as the estimated signals can be 
ordered differently and inversions are also possible. The overlapping sections of the 
estimated signals are ordered until this correlation matrix resembles a diagonal matrix. In 
this case, the signals will be in the same order. This technique is described in by G. 
Spence in "A Joint Estimation Approach for Periodic Signal Analysis", Queen's 
University Belfast, June 2003. Figure 14 illustrates this technique. In Figure 14(a) eight 
signals from one window are correlated with eight signals from an immediately following 
window. A peak 160 indicates that signal no. 8 from one window corresponds to signal 
no. 2 from the other window, and similar remarks apply to other peaks such as 162. 
Signal swapping is corrected by redesignating signals so that the same signal number is 
associated with correlated signals from the two windows. This yields the situation shown 
in Figure 14(b), which is as Figure 14(a) except for redesignation of signals. Here peaks 
such as 1 64 indicating correlated signals appear on the diagonal only. 
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Referring now to Figure 15, electrode apparatus suitable for implementing fetal ECG 
monitoring is indicated generally by 170. It comprises a number of electrodes 171a, 
171b, 171c, G and R suitable for placing on the surface of a mother's skin and 
monitoring vpltage signals developed there: here G indicates a ground electrode, R a 
5 reference electrode and 171a to 171b are signal electrodes. For convenience only five 
electrodes 171a etc.- are shown, but as many may be used as are necessary havi 
regard to the number of signals to be separated. The electrodes 171a etc. are 
connected via respective screened leads 172a to 172e to a lead box 173, which is also 
connected to a computer 1 74. 

10 The lead box 173 contains processing circuitry for signals from the signal electrodes 
171a etc. are shown inset at 173a: circuitry is shown for two signal electrodes 171a and 
171b, dots D indicate like circuitry for other signal electrodes. The circuitry comprises for 
each signal electrode 171a etc. a respective low-noise differential amplifier 175 with an 
analogue high-pass filter and an analogue low-pass anti-aliasing filter both incorporated 

15 in a filter unit 1 76. The filters 1 76 are connected to a common multi-channel analogue to 
digital converter (ADC) 177. The amplifiers 175 subtract a signal from the reference 
electrode R from signals from respective signal electrodes 171a etc., and amplify t 
resulting difference signals. The difference signals are converted to digital signals by 
ADO 177, and are then recorded and processed using the computer 174 to separate 

20 signals as described earlier. 
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Referring now to Figures 16 and 17, these show recording of signals and separated 
signals respectively for a further application of the invention, albeit physiological 
monitoring once more. The application is passive sleep recording to detect heart and 
breathing rates. This can help determine if a person is suffering from sleep apnoea 
25 (inability to breathe giving quickened heart rate). • 

Twelve signals 201 to 212 (marked on left) in Figure 16 were sampled at 250 Hz using 
12 sensors placed under a mattress on which a patient lay. The signals indicate vibration 
received by the sensocs though the bed/mattress. Fifty data windows each 4000 signal 
samples long (1 6 seconds) were processed using a window overtap of 85%. All available 
30 signals were processed, i.e. signals were not targeted. The reason for the large window 
overlap is that the invention is based on an assumption that conditions are slowly varying 
between successive windows, so the time difference between the beginnings of two 
successive windows must be sufficiently short to satisfy this assumption and make the 
initialisation process more accurate. In the present case, satisfying this criterion led to a 



31 



time difference of only 15% of a data window (2.4 seconds). When conditions pernnit. it is 
better to use non-overlapping contiguous windows as fewer data windows are required 
and processing tinne is nninimised. 

In Figure 17. twelve separated signals 221 to 232 (marked on left) are shown: df these 
signals 221. 224 and 225 are estimated heart signals with periodicity -0.83 seconds. 
Signals 226 and 228 are estimated signals associated with breathing with periodicity 
-4.4 seconds. 

It is possible in principle to envisage implementing signal separation without initialisation, 
and update processes, but with the above ordering technique to correct the problem of 
signal swapping. However, this has disadvantages. It would sacrifice the ability of the 
Invention to exploit the. computationally efficient targeted approach. Additionally, to 
provide reliable estimates in the cross-correlation technique, a large window overiap 
would have to be used. This means, that in order to process the data, a larger number of 
windows is required: e.g. with 50% window overlap the data is processed twice (and 
takes twice as long) compared to processing once (i.e. with no overlap). This is 
disadvantageous because it lengthens the time between detection of mixed signals by 
sensors. 14 and generation of unmixed signals. In the case of fetal ECG in particular, this 
is important because it may be necessary for a clinician to intervene rapidly if a heartbeat 
abnormality is detected, and it may be too late to save the fetus if processing time is too 
long. 

Finally, the invention can ajso be applied to complex-valued data. The principles of 
initialising the second and higher order stages and using small updates remain 
unchanged; however complex-valued versions of the second order stage (Jacobi) and 
higher order (ICA) stage are required. The use of ICA for complex-valued data is well 
researched and many decorrelation and ICA techniques have been developed for this 
problem. For more Information on a complex-valued Jacobi refer to: 

S. Haykin. "Adaptive filter theory, second edition", Prentice Hall, ISBN 0-13-013236-5, 
1991. 

In addition, more infomnation about complex-valued ICA methods can be found in: 

P. Comon, "Independent component analysis. A new concepf, Signal Processing 36, 
pp. 287-314. 1994. 
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IJ. Clarke. "Exploiting non-Gaussianity for signal separation". MILCOM 2000. 21st 
Century Military Communications Conference Proceedings. Volume: 2, 22-25 Oct. 2000. 



In this paper, Clarke showed that an extended version of BUSS can deal with complex 
signals by treating their- in-phase and quadrature components as individual real signals. 
By using Jacobi and this extended version of BLISS small updates of the initialis^ 
signals are again possible. This is because the smaller angle is chosen in the JacoDi 
technique and the rotation of the initialised signals in the higher order stage will again 

■ * 

depend on fourth order statistics. 

The equations for a complex-valued version of the invention, in the (k4-1)th window, can 
be summarised in four stages as follows. The superscript index H denotes the hemiitian 
of a matrix or vector. These equations can also be used to represent the real-valued 
case. 

1 . Initialisation in second order stage 



U*(X*.iX?,Ju, =U?(xJ,+A+iULi)Uit (39) 




2. Update using complex-valued Jacobi 




(40) 



(41) 



(42) 



3. Initialisation in iiiglier order stage 



(43) 



4. Update using complex-valued ICA 




(44) 
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Equations (41) and (42) represent updated eigenvectors and orthonormal signals 
respectively. Equations (44) and (45) represent updated independent signals and the 
higher rotation matrix respectively. The four stage process is summarised as described 
earlier with relerence to Figure 3. The targeted aspect of the invention is also applicable 
to the complex-valued case, and this is again involves updating only the signal pairs that 
include the desired signals. 

Equations given in the foregoing description for calculating intermediate quantities and 
results can clearly be evaluated by an appropriate computer program embodied m a 
carrier medium and running on a conventional computer system. Such a program is 
straightfonward for a skilled programmer to implement without requiring invention. Such a 
program and system will therefore not be described further. 
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Claims 



1. 



2. 



3. 



7. 



A method for dynamic blind signal- separation including processing signals 
associated with data windows characterised in that the method includes: ^ 

a) for pairs of successive data windows, each pair having a leading windSF 
and a following window, using results obtained in connection with each 
leading window to initialise data associated with the respective following 
window, and 

b) processing the initialised data to achieve signal separation. 

A method according to Claim 1 characterised in that it includes initialising both 
data orthogonality and data independence and doing so separately. 

A method according to Claim 1 characterised in that the data orthogonality is 
initialised and the initialised data is processed to update orthogonality using small 
updates to produce decorrelated data in a second order statistics procedure. 



4. A method according to Claim 3 characterised in that the decorrelated data is 
produced by a technique referred to as Jacobi and involving diagonalisation of a 
symmetric matrix by determining and applying rotations iteratively until off- 
diagonal elements of the matrix become substantially equal to zero. 

5. A method according to Clajm 3 characterised in that it includes a second stage of 
initialisation using results obtained for each leading window to initialise 
Independence of decorrelated data associated with the respective following 
window, this second stage using Independent component analysis (ICA) to apply 
small rotation updates to inifalised data in a higher than second order statistics 
procedure to produce signal independence and separation. 

6. • A method according to Claim 6 characterised in that the higher than second order 

statistics procedure is at least one of a third order and a fourth order statistics 

« • 

procedure. 




A method according to Claim 1 charafcterised in that it Is implemented in an 
acquisition phase in which signals are separated and desired signals are identified 
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among the separated signals, and in a subsequent phase in which only desired 
signals are processed to separation. 

Computer apparatus for dynamic blind signal separation programmed to process 
signals associated with data windows characterised in that the computer 
apparatus is also programmed to: 

a) process pairs of successive data windows, each pair having a leading 
window and a following window, 

b) use results obtained in connection with each leading window to initialise 
data associated with the respective following window, and 

c) process the initialised data to achieve signal separation. 

Computer apparatus according to Claim 8 characterised in that it is programmed 
to initialise both data orthogonality and data independence and to do so 
separately. 

Computer apparatus according to Claim 8 characterised in that it is programmed 
to initialise data orthogonality and process the initialised data to update 
orthogonality using small updates to produce decorrelated data in a second order 
statistics procedure. 

Computer apparatus according to Claim 10 characterised in that it is programmed 
to produce the decorrelated data by a technique referred to as JacobI and 
involving diagonalisation of a symmetric matrix by determining and applying 
rotations iteratlvely until off-diagonal elements of the matrix become substantially 
equal to zero. 

Computer apparatus according to Claim 10 characterised in that it is programmed 
to implement a second stage of initialisation using results obtained for each 
leading window to Initialise independence of decorrelated data associated with the 
respective following window, this second stage using ICA to apply small rotation 
updates to initialised data In a higher than second order statistics procedure to 
produce signal independence and separation. 
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Computer apparatus according to Claim 1 2 characterised in tliat the higher than 
second order statistics procedure is at least one of a third order and a fourth order 
statistics procedure. 

Computer apparatus according to Claim 8 characterised in that it is programm^^ 
to implement an acquisition phase in which signals are separated and desired 
signals are identified among the separated signals, and a subsequent phase in 
which only desired signals are processed to separation. 

Computer software for dynamic blind signal separation including processing 
signals associated with data windows characterised in that the software has 
instructions for controlling computer apparatus to: 

a) process pairs of successive data windows, each pair having a. leading 
■ window and a following window, 

b) use results obtained in connection with each leading window to initialise 
data associated with the respective following window, and 

c) process the initialised data to achieve signal separation. 




Computer software according to Claim 15 characterised in that it includes 
instructions for implementing initialisation of both data orthogonality and data 
independence and for doing so separately. 

Computer software according to Claim 15 characterised in that it includes 
instructions for initialising the data orthogonality and for processing the initialised 
data to update orthogonality using small updates to produce decorrelated data in 
a second order statistics procedure. 

Computer software according to Claim 17 characterised in that it includes 
instructions for producing the decorrelated data by a technique referred to as 
Jacobi and involving diagonalisation of a symmetric matrix by determining and 
applying rotations iteratively until off-diagonal elements of the matrix become 
substantially equal to zero. 

Computer software according to Claim 17 characterised in that it includes 
instructions for implementing a second stage of initialisation using results obtained 
for each leading window to initialise independence of decorrelated data 
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associated with tiie respective following window, this second stage using 
independent component analysis (ICA) to apply small rotation updates to 
initialised data in a higher than second order statistics procedure to produce signal 
independence and separation. 

Computer software according to Claim 19 characterised in that the higher than 
second order statistics procedure is at least one of a third order and a fourth order 
statistics procedure. 

Computer software according to Claim 1 5 characterised in that it is implemented 
in an acquisition phase in which signals are separated' and desired signals are 
identified among the separated signals, and in a subsequent phase in which only 
desired signals are processed to separation. 
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Abstract 



A method for dynamic blind signal separation generates initialisation information by 
processing an immediately preceding data window, This information is input at 21 and 
used at 22 to initialise orthogonality of data in an immediately following window. Initialis^ 
data are decorrelated at 23 with small update angles using a Jacobi technique. Steps 22 
and 23 are collectively a second order stage of processing in statistical terms producing 
orthonormal signals. The orthonormal signals are initialised at 25 and then undergo 
separation at 27 by ICA with small angle updates using higher order statistics, i.e. higher 
than second order. Steps 25 and 27 are collectively a higher order stage of processing in 
statistical terms producing separated signals at an output 28. The method may be 
implemented in an acquisition phase to separate signals and among them identify desired 
signals, and a subsequent phase in which only the desired signals are separated. 



Figure 3 should accompany the Abstract. 
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